/* START 2_tripCountsByDayInPanel.do */

/* this script compares trip counts by day for the panel data by
origin/destination for the descriptive section. The goal is to compare the raw
number of trips before/after the event break. We will use JUN 2, 2015 as a
benchmark (it is a Monday) */

/* prerequisite packages: fillmissing, labutil2, gtools, parallel, sepscatter */

version 16.1
set more off
cap clear all

/* PATHS */
use "code/config.do"

/* Initialize 12 threads */
// parallel initialize 12

/* load accessory functions */
qui do "./code/99_accessoryFunctions.do"
use "data/panel.dta", clear

/* fix labeling from earlier */
cap label copy oa1 oa
cap label drop oa1
cap label values oa da oa

/*****************************/
/** COLLAPSE THE PANEL DATA **/
/*****************************/

/* compute trip dates */
gen tripDate = dofc(START)
format tripDate %td
drop if tripDate == .
gen byte after = (tripDate >= td(27dec2015))

/* finer time subdivisions */
gen byte hr = hh(START)

/* destinations */
gen byte oCCK = (oa == "JELEBU":oa) | (oa == "HILLVIEW":oa) | (oa == "GOMBAK":oa)
gen byte dCCK = (da == "JELEBU":oa) | (da == "HILLVIEW":oa) | (da == "GOMBAK":oa)
gegen byte gCCK = group(oCCK dCCK)
gen byte CCK = (oCCK == 1) | (dCCK == 1)
gen byte obw = (oa == "ANAK BUKIT":oa)
gen byte dbw = (da == "ANAK BUKIT":oa)
gen byte bw = (obw == 1) | (dbw == 1)
gegen byte gbw = group(obw dbw)

gen byte odt = (op == "DOWNTOWN CORE":op)
gen byte ddt = (dp == "DOWNTOWN CORE":op)
gen byte dt = (odt == 1) | (ddt == 1)
gegen byte gdt = group(odt ddt)

gen byte geylang = (op == "GEYLANG":op) | (dp == "GEYLANG":op)
gen byte nearDTL = (oa == "GOMBAK":oa | da == "GOMBAK":oa) |  ///
  (oa == "JELEBU":oa | da == "JELEBU":oa) |  ///
  (oa == "HILLVIEW":oa | da == "HILLVIEW":oa) |  ///
  (oa == "ANAK BUKIT":oa | da == "ANAK BUKIT":oa) |  ///
  (oa == "HOLLAND ROAD":oa | da == "HOLLAND ROAD":oa) |  ///
  (oa == "SWISS CLUB":oa | da == "SWISS CLUB":oa) |  ///
  (oa == "CORONATION ROAD":oa | da == "CORONATION ROAD":oa) |  ///
  (oa == "TYERSALL":oa | da == "TYERSALL":oa) |  ///
  (oa == "NASSIM":oa | da == "NASSIM":oa) |  ///
  (oa == "NEWTON CIRCUS":oa | da == "NEWTON CIRCUS":oa) |  ///
  (oa == "FARRER PARK":oa | da == "FARRER PARK":oa) |  ///
  (oa == "BENCOOLEN":oa | da == "BENCOOLEN":oa) |  ///
  (oa == "BUGIS":oa | da == "BUGIS":oa)

/* TRIP COUNTS: CCK vs. NOT CCK */
frames put oa CCK gCCK tripDate hr TKT_TYP_CD, into(lrcck)
fv lrcck
preserve
gcollapse (count) oa, by(CCK tripDate)
ren oa trips
gen lastday = (day(tripDate) == 8)
drop if lastday == 1

gegen tripRef = mean(trips) if tripDate == td(02jun2015), by(CCK)
fillmissing tripRef, by(CCK)
gen normTrips = trips/tripRef

binscatter2 normTrips tripDate, ///
  by(CCK) linetype(connect) xla(, format(%td)) ///
  xline(20449, lcolor(green) lpattern(dash)) ///
  xlabel(#4, labsize(vsmall)) ylabel(, labsize(vsmall)) ///
  title("Binned scaled average daily trip counts" ///
  "to/from CCK or not") ///
  xtitle("Date") ytitle("Avg. trip count (rel. Jun 2 2015)") ///
  legend(order(1 "All other trips (Jun 2: 4.34m)" ///
  2 "To/From/Within Gombak, Hillview, Jelebu (Jun 2: 57.7k)") rows(2))
graph export "${make_data}/aggNormTripsToFromCCKPanel.png", replace width(1920)

restore

/* TRIP COUNTS BY BY CARD TYPE */
preserve
gcollapse (count) oa, by(CCK tripDate TKT_TYP_CD)
ren oa trips
gen lastday = (day(tripDate) == 8)
drop if lastday == 1

gegen grp = group(TKT_TYP_CD CCK)
gen logTrips = log(trips)
gegen tripRef = mean(trips) if tripDate == td(02jun2015), by(grp)
fillmissing tripRef, by(grp)
gen normTrips = trips/tripRef

/* raw numbers */
binscatter2 normTrips tripDate if TKT_TYP_CD == 34 | TKT_TYP_CD == 38, ///
  by(grp) linetype(connect) xla(, format(%td)) ///
  xline(20449, lcolor(purple) lpattern(dash)) ///
  xlabel(#4, labsize(vsmall)) ylabel(, labsize(vsmall)) ///
  title("Binned scaled average daily trip counts by card type" ///
  "to/from CCK or not") ///
  xtitle("Date") ytitle("Avg. trip count (rel. Jun 2 2015)") ///
  legend(order(1 "Adult farecard trips outside CCK (Jun 2: 2.73m)" ///
  2 "Adult farecard trips to/from/within CCK (Jun 2: 32.3k)" ///
  3 "Subsidized worker trips outside CCK (Jun 2: 212.9k)" ///
  4 "Subsidized worker trips to/from/within CCK (Jun 2: 3.16k)") rows(4))
/* graph export "${make_data}/norm3438TripsCCKPanel.png", replace width(1920) */
graph save "${dropbox}/norm3438TripsCCKPanel.gph", replace

binscatter2 normTrips tripDate if TKT_TYP_CD >= 40 & TKT_TYP_CD <= 42, ///
  by(grp) linetype(connect)  ///
  msymbols(o oh d dh t th) ///
  xline(20449, lcolor(purple) lpattern(dash)) ///
  xla(, format(%td) labsize(vsmall)) ylabel(, labsize(vsmall)) ///
  title("Binned scaled avg. daily trip counts " ///
  "by student card type & to/from CCK or not") ///
  xtitle("Date") ytitle("Avg. trip count (rel. Jun 2 2015)") ///
  legend(order(1 "Elem., not (Jun 2: 87.8k)" ///
  2 "Elem., CCK (Jun 2: 1.52k)" ///
  3 "Mid., not (Jun 2: 243.1k)" ///
  4 "Mid., CCK (Jun 2: 4.68k)" ///
  5 "Tertiary, not (Jun 2: 60.8k)" ///
  6 "Tertiary, CCK (Jun 2: 860)") rows(3))
graph export "${make_data}/norm40-42TripsCCKPanel.png", replace width(1920)

// trip counts by location
keep normTrips TKT_TYP_CD tripDate CCK
greshape wide normTrips, by(TKT_TYP_CD tripDate) keys(CCK)
gen dNormTrips = normTrips1 - normTrips0

binscatter2 dNormTrips tripDate if TKT_TYP_CD == 34 | TKT_TYP_CD == 38, ///
  by(TKT_TYP_CD) linetype(connect) xla(, format(%td)) ///
  xline(20449, lcolor(purple) lpattern(dash)) ///
  xlabel(#4, labsize(vsmall)) ylabel(, labsize(vsmall)) ///
  msymbols(o oh d dh) ///
  title("Differences in binned scaled average daily trip counts," ///
  "to/from/within CCK vs. not, by card type") ///
  xtitle("Date") ytitle("Differences in average trip count (base: Jun 2 2015)") ///
  legend(order(1 "Adult farecard trips" ///
  2 "Subsidized worker trips") rows(2))
graph save "${dropbox}/normDiff3438TripsCCKPanel.gph", replace
restore

/* TRIP COUNTS BY HOUR, BEFORE VS AFTER */
preserve
gcollapse (count) n = oa, by(gCCK tripDate hr)
gen byte after = (tripDate >= td(27dec2015))
gen logn = log(n)
gegen grp = group(gCCK after)
binscatter2 logn hr, by(grp) linetype(none)

gegen nRef = mean(n) if tripDate == td(02jun2015), by(gCCK)
fillmissing nRef, by(gCCK)
gen nnorm = n/nRef
/* gen lognnorm = log(nnorm) */
binscatter2 nnorm hr, by(grp) linetype(connect) ///
  msymbols(o oh d dh t th s sh) ///
  title("Scaled trip counts by hour of day by origin & destination," ///
  "in CCK (GHJ), 06/2015-06/2018") ///
  xlabel(#12, labsize(small)) ylabel(, labsize(small)) ///
  xtitle("Hour") ///
  ytitle("Trip count (rel. grp. mean on Jun 2 2015)") ///
  legend(order(1 "Outside CCK, pre-DTL2" ///
  2 "Outside CCK, post-DTL2" ///  
  3 "Starting in CCK, pre-DTL2" ///
  4 "Starting in CCK, post-DTL2" ///
  5 "Ending in CCK, pre-DTL2" ///
  6 "Ending in CCK, post-DTL2" ///
  7 "Within CCK, pre-DTL2" ///
  8 "Within CCK, post-DTL2" ///
  ) rows(4))
graph export "${make_data}/aggNormODTripsPerHourCCKPanel.png", replace width(1920)
restore

// repeat by card type and direction of travel
preserve
gcollapse (count) n = oa, by(gCCK tripDate hr TKT_TYP_CD)
gen byte after = (tripDate >= td(27dec2015))
/* gen logn = log(n) */
gegen grp = group(gCCK TKT_TYP_CD after)
gegen nRef = mean(n) if tripDate == td(02jun2015), by(gCCK TKT_TYP_CD)
fillmissing nRef, by(gCCK TKT_TYP_CD)
gen nnorm = n/nRef
/* gen lognnorm = log(nnorm) */
foreach ttype in 34 38 39 40 41 42 {
	qui binscatter nnorm hr if TKT_TYP_CD == `ttype' & hr >= 5, by(grp) linetype(connect) ///
	  msymbols(o oh d dh t th s sh) ///
	  title("Scaled trip counts by hour of day by origin & destination," ///
	  "in CCK (GHJ), 06/2015-06/2018 (Farecard ID: `ttype')") ///
	  xlabel(#12, labsize(small)) ylabel(, labsize(small)) ///
	  xtitle("Hour") ///
	  ytitle("Trip count (rel. grp. mean on Jun 2 2015)") ///
	  legend(order(1 "Outside CCK, pre-DTL2" ///
	  2 "Outside CCK, post-DTL2" ///  
	  3 "Starting in CCK, pre-DTL2" ///
	  4 "Starting in CCK, post-DTL2" ///
	  5 "Ending in CCK, pre-DTL2" ///
	  6 "Ending in CCK, post-DTL2" ///
	  7 "Within CCK, pre-DTL2" ///
	  8 "Within CCK, post-DTL2" ///
	  ) rows(4))
	qui graph export "${make_data}/norm`ttype'ODTripsPerHourCCKPanel.png", replace width(1920)
	}

gcollapse (mean) n, by(hr gCCK TKT_TYP_CD after)
greshape wide n, by(gCCK TKT_TYP_CD hr) keys(after)
gen ndiff = n1-n0
gen fracdiff = ndiff/n0
gegen grp = group(gCCK TKT_TYP_CD)
foreach ttype in 34 38 39 40 41 42 {
	qui sepscatter fracdiff hr ///
	  if (TKT_TYP_CD == `ttype') & hr >= 5, ///
	  separate(grp) recast(connect) ///
	  lpattern(solid solid shortdash shortdash) ///
	  xlabel(#12, labsize(small)) ylabel(, labsize(small)) ///
	  title("Scaled trip count diffs by orig./dest. & time of day," ///
	  "to CCK (GHJ), 06/2015-06/2018 (Farecard ID: `ttype')") ///
	  xtitle("Hour") ytitle("Trip count diff. (rel. pre-DTL2)") ///
	  legend(order(1 "Outside CCK" ///
	  2 "Starting in CCK" ///  
	  3 "Ending in CCK" ///
	  4 "Within CCK" ///
	  ) rows(2))
	qui graph export "${make_data}/normDiff`ttype'ODTripsPerHourCCKPanel.png", replace width(1920)
	}
restore
fv

/* BENCHMARKS FOR BEAUTY WORLD */

frames put oa bw gbw tripDate hr TKT_TYP_CD, into(lrbw)
fv lrbw
preserve
gcollapse (count) trips = oa, by(bw tripDate TKT_TYP_CD)
gen lastday = (day(tripDate) == 8)
drop if lastday == 1

gegen grp = group(TKT_TYP_CD bw)
gen logTrips = log(trips)
gegen tripRef = mean(trips) if tripDate == td(02jun2015), by(grp)
fillmissing tripRef, by(grp)
gen normTrips = trips/tripRef

binscatter2 normTrips tripDate if TKT_TYP_CD == 34 | TKT_TYP_CD == 38, ///
  by(grp) linetype(connect) xla(, format(%td)) ///
  xline(20449, lcolor(purple) lpattern(dash)) ///
  xlabel(#4, labsize(vsmall)) ylabel(, labsize(vsmall)) ///
  msymbols(o oh d dh) ///
  title("Binned scaled average daily trip counts by card type" ///
  "to/from Beauty World or not") ///
  xtitle("Date") ytitle("Avg. trip count (rel. Jun 2 2015)") ///
  legend(order(1 "Adult farecard trips outside BW (Jun 2: 2.74m)" ///
  2 "Adult farecard trips to/from/within BW (Jun 2: 18.7k)" ///
  3 "Subsidized worker trips outside BW (Jun 2: 215.1k)" ///
  4 "Subsidized worker trips to/from/within BW (Jun 2: 1.04k)") rows(4))
/* graph export "${make_data}/norm3438TripsBWPanel.png", replace width(1920) */
graph save "${dropbox}/norm3438TripsBWPanel.gph", replace


binscatter2 normTrips tripDate if TKT_TYP_CD >= 40 & TKT_TYP_CD <= 42, ///
  by(grp) linetype(connect)  ///
  msymbols(o oh d dh t th) ///
  xline(20449, lcolor(purple) lpattern(dash)) ///
  xla(, format(%td) labsize(vsmall)) ylabel(, labsize(vsmall)) ///
  title("Binned scaled avg. daily trip counts " ///
  "by student card type & to/from BW or not") ///
  xtitle("Date") ytitle("Avg. trip count (rel. Jun 2 2015)") ///
  legend(order(1 "Elem., not (Jun 2: 88.7k)" ///
  2 "Elem., BW (Jun 2: 631)" ///
  3 "Mid., not (Jun 2: 245.8k)" ///
  4 "Mid., BW (Jun 2: 1.93k)" ///
  5 "Tertiary, not (Jun 2: 61.2k)" ///
  6 "Tertiary, BW (Jun 2: 446)") rows(3))
graph export "${make_data}/norm40-42TripsBWPanel.png", replace width(1920)
restore

// beauty world, by hour
preserve
gcollapse (count) n = oa, by(gbw tripDate hr)
gen byte after = (tripDate >= td(27dec2015))
/* gen logn = log(n) */
gegen grp = group(gbw after)
/* binscatter2 logn hr, by(grp) linetype(none) */

gegen nRef = mean(n) if tripDate == td(02jun2015), by(gbw)
fillmissing nRef, by(gbw)
gen nnorm = n/nRef
/* gen lognnorm = log(nnorm) */
binscatter nnorm hr, by(grp) linetype(connect) ///
  msymbols(o oh d dh t th s sh) ///
  title("Scaled trip counts by hour of day by origin & destination," ///
  "in Beauty World, 06/2015-06/2018") ///
  xlabel(#12, labsize(small)) ylabel(, labsize(small)) ///
  xtitle("Hour") ///
  ytitle("Trip count (rel. grp. mean on Jun 2 2015)") ///
  legend(order(1 "Outside BW, pre-DTL2" ///
  2 "Outside BW, post-DTL2" ///  
  3 "Starting in BW, pre-DTL2" ///
  4 "Starting in BW, post-DTL2" ///
  5 "Ending in BW, pre-DTL2" ///
  6 "Ending in BW, post-DTL2" ///
  7 "Within BW, pre-DTL2" ///
  8 "Within BW, post-DTL2" ///
  ) rows(4))
graph export "${make_data}/aggNormODTripsPerHourBWPanel.png", replace width(1920)
restore

// beauty world, by card type and hour of day
preserve
gcollapse (count) n = oa, by(gbw tripDate hr TKT_TYP_CD)
gen byte after = (tripDate >= td(27dec2015))
/* gen logn = log(n) */
gegen grp = group(gbw TKT_TYP_CD after)
gegen nRef = mean(n) if tripDate == td(02jun2015), by(gbw TKT_TYP_CD)
fillmissing nRef, by(gbw TKT_TYP_CD)
gen nnorm = n/nRef
/* gen lognnorm = log(nnorm) */
foreach ttype in 34 38 39 40 41 42 {
	qui binscatter nnorm hr if TKT_TYP_CD == `ttype' & hr >= 5, by(grp) linetype(connect) ///
	  msymbols(o oh d dh t th s sh) ///
	  title("Scaled trip counts by hour of day by origin & destination," ///
	  "in Beauty World, 06/2015-06/2018 (Farecard ID: `ttype')") ///
	  xlabel(#12, labsize(small)) ylabel(, labsize(small)) ///
	  xtitle("Hour") ///
	  ytitle("Trip count (rel. grp. mean on Jun 2 2015)") ///
	  legend(order(1 "Outside BW, pre-DTL2" ///
	  2 "Outside BW, post-DTL2" ///  
	  3 "Starting in BW, pre-DTL2" ///
	  4 "Starting in BW, post-DTL2" ///
	  5 "Ending in BW, pre-DTL2" ///
	  6 "Ending in BW, post-DTL2" ///
	  7 "Within BW, pre-DTL2" ///
	  8 "Within BW, post-DTL2" ///
	  ) rows(4))
	qui graph export "${make_data}/norm`ttype'ODTripsPerHourBWPanel.png", replace width(1920)
	}

gcollapse (mean) n, by(hr gbw TKT_TYP_CD after)
greshape wide n, by(gbw TKT_TYP_CD hr) keys(after)
gen ndiff = n1-n0
gen fracdiff = ndiff/n0
gegen grp = group(gbw TKT_TYP_CD)
foreach ttype in 34 38 39 40 41 42 {
	qui sepscatter fracdiff hr ///
	  if (TKT_TYP_CD == `ttype') & hr >= 5, ///
	  separate(grp) recast(connect) ///
	  lpattern(solid solid shortdash shortdash) ///
	  xlabel(#12, labsize(small)) ylabel(, labsize(small)) ///
	  title("Scaled trip count diffs by orig./dest. & time of day," ///
	  "to Beauty World, 06/2015-06/2018 (Farecard ID: `ttype')") ///
	  xtitle("Hour") ytitle("Trip count diff. (rel. pre-DTL2)") ///
	  legend(order(1 "Outside BW" ///
	  2 "Starting in BW" ///  
	  3 "Ending in BW" ///
	  4 "Within BW" ///
	  ) rows(2))
	qui graph export "${make_data}/normDiff`ttype'ODTripsPerHourBWPanel.png", replace width(1920)
	}
restore
fv

/* BENCHMARKS FOR DOWNTOWN */
frames put oa dt gdt tripDate hr TKT_TYP_CD, into(lrdt)
fv lrdt

// downtown, by hour
preserve
gcollapse (count) n = oa, by(gdt tripDate hr)
gen byte after = (tripDate >= td(27dec2015))
/* gen logn = log(n) */
gegen grp = group(gdt after)
/* binscatter2 logn hr, by(grp) linetype(none) */

gegen nRef = mean(n) if tripDate == td(02jun2015), by(gdt)
fillmissing nRef, by(gdt)
gen nnorm = n/nRef
/* gen lognnorm = log(nnorm) */
binscatter nnorm hr if hr >= 5, by(grp) linetype(connect) ///
  msymbols(o oh d dh t th s sh) ///
  title("Scaled trip counts by hour of day by origin & destination," ///
  "Downtown, 06/2015-06/2018") ///
  xlabel(#12, labsize(small)) ylabel(, labsize(small)) ///
  xtitle("Hour") ///
  ytitle("Trip count (rel. grp. mean on Jun 2 2015)") ///
  legend(order(1 "Outside DT, pre-DTL2" ///
  2 "Outside DT, post-DTL2" ///  
  3 "Starting in DT, pre-DTL2" ///
  4 "Starting in DT, post-DTL2" ///
  5 "Ending in DT, pre-DTL2" ///
  6 "Ending in DT, post-DTL2" ///
  7 "Within DT, pre-DTL2" ///
  8 "Within DT, post-DTL2" ///
  ) rows(4))
graph export "${make_data}/aggNormODTripsPerHourDTPanel.png", replace width(1920)
restore

// downtown, by card type and hour of day
preserve
gcollapse (count) n = oa, by(gdt tripDate hr TKT_TYP_CD)
gen byte after = (tripDate >= td(27dec2015))
/* gen logn = log(n) */
gegen grp = group(gdt TKT_TYP_CD after)
gegen nRef = mean(n) if tripDate == td(02jun2015), by(gdt TKT_TYP_CD)
fillmissing nRef, by(gdt TKT_TYP_CD)
gen nnorm = n/nRef
/* gen lognnorm = log(nnorm) */
foreach ttype in 34 38 39 40 41 42 {
	qui binscatter nnorm hr if TKT_TYP_CD == `ttype' & hr >= 5, by(grp) linetype(connect) ///
	  msymbols(o oh d dh t th s sh) ///
	  title("Scaled trip counts by hour of day by origin & destination," ///
	  "Downtown, 06/2015-06/2018 (Farecard ID: `ttype')") ///
	  xlabel(#12, labsize(small)) ylabel(, labsize(small)) ///
	  xtitle("Hour") ///
	  ytitle("Trip count (rel. grp. mean on Jun 2 2015)") ///
	  legend(order(1 "Outside DT, pre-DTL2" ///
	  2 "Outside DT, post-DTL2" ///  
	  3 "Starting in DT, pre-DTL2" ///
	  4 "Starting in DT, post-DTL2" ///
	  5 "Ending in DT, pre-DTL2" ///
	  6 "Ending in DT, post-DTL2" ///
	  7 "Within DT, pre-DTL2" ///
	  8 "Within DT, post-DTL2" ///
	  ) rows(4))
	qui graph export "${make_data}/norm`ttype'ODTripsPerHourDTPanel.png", replace width(1920)
	}

gcollapse (mean) n, by(hr gdt TKT_TYP_CD after)
greshape wide n, by(gdt TKT_TYP_CD hr) keys(after)
gen ndiff = n1-n0
gen fracdiff = ndiff/n0
gegen grp = group(gdt TKT_TYP_CD)
foreach ttype in 34 38 39 40 41 42 {
	qui sepscatter fracdiff hr ///
	  if (TKT_TYP_CD == `ttype') & hr >= 5, ///
	  separate(grp) recast(connect) ///
	  lpattern(solid solid shortdash shortdash) ///
	  xlabel(#12, labsize(small)) ylabel(, labsize(small)) ///
	  title("Scaled trip count diffs by orig./dest. & time of day," ///
	  "Downtown, 06/2015-06/2018 (Farecard ID: `ttype')") ///
	  xtitle("Hour") ytitle("Trip count diff. (rel. pre-DTL2)") ///
	  legend(order(1 "Outside DT" ///
	  2 "Starting in DT" ///  
	  3 "Ending in DT" ///
	  4 "Within DT" ///
	  ) rows(2))
	qui graph export "${make_data}/normDiff`ttype'ODTripsPerHourDTPanel.png", replace width(1920)
	}
restore
fv

/* CONTROL: GEYLANG */
preserve
gcollapse (count) oa, by(geylang tripDate TKT_TYP_CD)
ren oa trips
gen lastday = (day(tripDate) == 8)
drop if lastday == 1

gegen grp = group(TKT_TYP_CD geylang)
gen logTrips = log(trips)
gegen tripRef = mean(trips) if tripDate == td(02jun2015), by(grp)
fillmissing tripRef, by(grp)
gen normTrips = trips/tripRef

binscatter2 normTrips tripDate if TKT_TYP_CD == 34 | TKT_TYP_CD == 38, ///
  by(grp) linetype(connect) xla(, format(%td)) ///
  xline(20449, lcolor(purple) lpattern(dash)) ///
  xlabel(#4, labsize(vsmall)) ylabel(, labsize(vsmall)) ///
  msymbols(o oh d dh) ///
  title("Binned scaled average daily trip counts by card type" ///
  "to/from Geylang or not") ///
  xtitle("Date") ytitle("Avg. trip count (rel. Jun 2 2015)") ///
  legend(order(1 "Adult farecard trips outside Geylang (Jun 2: 2.52m)" ///
  2 "Adult farecard trips to/from/within Geylang (Jun 2: 241.8k)" ///
  3 "Subsidized worker trips outside Geylang (Jun 2: 197.2k)" ///
  4 "Subsidized worker trips to/from/within Geylang (Jun 2: 18.9k)") rows(4))
graph export "${make_data}/norm3438TripsGeylangPanel.png", replace width(1920)

binscatter2 normTrips tripDate if TKT_TYP_CD >= 40 & TKT_TYP_CD <= 42, ///
  by(grp) linetype(connect)  ///
  msymbols(o oh d dh t th) ///
  xline(20449, lcolor(purple) lpattern(dash)) ///
  xla(, format(%td) labsize(vsmall)) ylabel(, labsize(vsmall)) ///
  title("Binned scaled avg. daily trip counts " ///
  "by student card type & to/from Geylang or not") ///
  xtitle("Date") ytitle("Avg. trip count (rel. Jun 2 2015)") ///
  legend(order(1 "Elem., not (Jun 2: 85.3k)" ///
  2 "Elem., Geyl. (Jun 2: 4.01k)" ///
  3 "Mid., not (Jun 2: 236.5k)" ///
  4 "Mid., Geyl. (Jun 2: 11.3k)" ///
  5 "Tertiary, not (Jun 2: 58.0k)" ///
  6 "Tertiary, Geyl. (Jun 2: 3.70k)") rows(3))
graph export "${make_data}/norm40-42TripsGeylangPanel.png", replace width(1920)


/* LR DTL2 binscatter (diffs) */
preserve
gcollapse (count) trips = oa, by(nearDTL tripDate TKT_TYP_CD)
gen lastday = (day(tripDate) == 8)
drop if lastday == 1

gegen grp = group(TKT_TYP_CD nearDTL)
gen logTrips = log(trips)
gegen tripRef = mean(trips) if tripDate == td(02jun2015), by(grp)
fillmissing tripRef, by(grp)
gen normTrips = trips/tripRef

/* binscatter2 normTrips tripDate if TKT_TYP_CD == 34 | TKT_TYP_CD == 38, /// */
/*   by(grp) linetype(connect) xla(, format(%td)) /// */
/*   xline(20449, lcolor(purple) lpattern(dash)) /// */
/*   xlabel(#4, labsize(vsmall)) ylabel(, labsize(vsmall)) /// */
/*   msymbols(o oh d dh) /// */
/*   title("Binned scaled average daily trip counts by card type" /// */
/*   "to/from Beauty World or not") /// */
/*   xtitle("Date") ytitle("Avg. trip count (rel. Jun 2 2015)") /// */
/*   legend(order(1 "Adult farecard trips outside BW (Jun 2: 2.74m)" /// */
/*   2 "Adult farecard trips to/from/within BW (Jun 2: 18.7k)" /// */
/*   3 "Subsidized worker trips outside BW (Jun 2: 215.1k)" /// */
/*   4 "Subsidized worker trips to/from/within BW (Jun 2: 1.04k)") rows(4)) */

keep normTrips TKT_TYP_CD tripDate nearDTL
greshape wide normTrips, by(TKT_TYP_CD tripDate) keys(nearDTL)
gen dNormTrips = normTrips1 - normTrips0

binscatter2 dNormTrips tripDate if TKT_TYP_CD == 34 | TKT_TYP_CD == 38, ///
  by(TKT_TYP_CD) linetype(connect) xla(, format(%td)) ///
  xline(20449, lcolor(purple) lpattern(dash)) ///
  xlabel(#4, labsize(vsmall)) ylabel(, labsize(vsmall)) ///
  msymbols(o oh d dh) ///
  title("Differences in binned scaled average daily trip counts," ///
  "to/from/within DTL2 subzones vs. not, by card type") ///
  xtitle("Date") ytitle("Differences in average trip count (base: Jun 2 2015)") ///
  legend(order(1 "Adult farecard trips" ///
  2 "Subsidized worker trips") rows(2))
graph save "${dropbox}/normDiff3438TripsDTL2Panel.gph", replace
